Why?

Deep Dive into quantile residuals for model evaluation and the most popular implementation of them, DHARMa.

State of DHARMa

• DHARMa was first released in 2016. CRAN downloads were sparse until 2019, and began steadily increasing. It averages approximately ~500 downloads per day (source: https://hadley.shinyapps.io/cran-downloads/)

• It has reverse dependencies with 18 R packaged including glmmTMB, performance, easystats

• DHARMa does not have a published peer reviewed manuscript paper supporing its release

Standard Residuals from GLMMs Can Be Nonsensical

Linear Mixed Model

\[Y_{ij} = X_i\beta + Z_ja + \epsilon_{ij} \]

\[\epsilon \sim N(0, \sigma^2 \mathbf{I_n}) \] \[a \sim N(0, \sigma_a^2) \]

\[ y_i|a_{j} \sim 𝑁(x_i \beta + a_j, \sigma^2) \]

\[ y_i \sim 𝑁(x_i \beta, \sigma^2) \]

‘Standard’ Residuals

Raw \[\epsilon_i = Y_i - \hat{Y_i} \]

(internally) Studentized:

\[\epsilon_i = \frac {Y_i - \hat{Y_i}} {\hat{s}}\]

Pearson/Scaled:

\[ \epsilon_i = \frac {Y_i - \hat{Y_i}} {sd(Y_i)} \]

GLMM Formulation: Binomial Example

Dependent variable: \(Y_{ij}\) (e.g. survival)

survival: \(p_i = Y_i/N\)

\[Y_{ij}|r_j \sim Binomial(N, \pi_{ij})\]

\[\eta_{ij} = \eta + \alpha_i + r_j\] (\(\eta\) = mean, \(\alpha_i\) = treatment, \(r_j\) = grouping variable)

\[\eta_{ij} = log(\frac{\pi_{ij}} {1 - \pi_{ij}} )\]

Raw/Scaled Residuals Can Lack Diagnostic Capabilities for GLMMs

• Visual patterns are difficult to interpret
• Over and under-dispersion are difficult to assess
• Goodness of fit tests are not valid outside of normal-distriubuted variables

Different distributional assumptions and mathematical mean/variance relationships of some distributions require a different approach

Background

For \(y_1,...,y_n\) responses:

\[y_i \sim \mathcal{P}(\mu_i, \phi)\]

CDF: \[F(y; \mu, \phi)\]

For a continuous \(F\), \(F(y_i; \mu_i, \phi)\) are uniformly distributed, and the quantile residuals are:

\[ r_{i,q} = \Phi^{-1} \left\{F(y_i; \hat{\mu_i}, \hat{\phi} )\right\}\]

Example: Survival Times of Patients

Survival: \(y_i \sim Exp(\mu_i)\)

\[ \text{log }\mu_i = \beta_0 + \beta_1 \text{ log } x_i \]

Quantile residuals:

\[ r_{i} = \Phi^{-1} \left\{ 1-exp(y_i/\hat{\mu_i}) \right\}\]

Discontinuous CDF, F

\[r_i = \Phi^{-1}(u_i)\] \[ u_i \sim Uniform(a_i, b_j]\] \[a_i = lim_{y \rightarrow y_i} F(y; \hat{\mu_i}, \hat{\phi})\]

\[ b_i = F(y_i; \hat{\mu_i}, \hat{\phi}) \]

Method implemented in ‘statmod’

This Method Has Been Around

The DHARMa Process

1. Model something using a generalized linear model with a chosen distribution and link function.

\[\mathbf{Y = X\beta+Za}\] \[ \mathbf{Y|a\sim G(\mu, R)} \] linear predictor: \(\mathbf{\eta = X\beta}\)

fitted value: \(E(Y) = \mu\)

Link function: \(\eta = g(\mu)\)

The DHARMa Process

For each observation in the data set:

2. Simulate new observations (\(n_{sim} = n_{data}\)) using fitted values as the model distributional parameters (e.g. shape, scale) as appropriate. \[ Y_i \sim G(\hat{\mu_i}, \hat{\phi}) \]

3. Calculate the quantile for the cumulative distribution function

\[ r_{i} = F(y_i; \hat{\mu_i}, \hat{\phi} )\]

Expectations of the Quantile Residuals

\[ r_{ij} \sim Uniform(0,1) \]

\(r_{ij} = 0\): everything is larger

\(r_{ij} = 1\): everything is smaller

\(r_{ij} = 0.5\): right in the middle

DHARMa runs 250 simulations (per observation) by default, they recommend up to 1000

Poisson Example

\[ (\hat{Y_i} = \lambda = 5,; \quad Y_i = 7 )\]

Quantile Residuals for 1 Observation

Poisson GLMM Example

Create a data set with a random effect (“group”, 10 levels), a covariate (“Environment1”) and a Poisson-distributed response variable (“observedResponse”):

testData = createData(sampleSize = 500)

Analyze correctly and incorrectly:

rightModel <- glmer(observedResponse ~ Environment1 + (1|group) , 
                     family = "poisson", data = testData)

wrongModel <- lmer(observedResponse ~ Environment1 + (1|group) , 
                     data = testData)

‘Standard’ Residuals

Quantile Residuals

Correctly Specified Model

Quantile Residuals

Incorrectly Specified Model

Distributions of Quantile Residuals

Quantile Function of the Standard Normal Distribution

Notes on DHARMa Implementation

• It depends on the simulation functions build into GLMM package
• Supported packages: lm(), `glm(), lme4, mgcv, glmmTMB spaMM, GLMMadaptive, brms, & more
• Commonly, only the last stochastic level (e.g. Poisson) is simulated, conditional on the fitted random effects – basically an omnibus test
• Residuals can be simulated for individuals predictors and tests can be conducted for individual factors

lme4::glmer() Marginal Model

By default, simulations are conducted on a marginal model:

rightModel <- glmer(observedResponse ~ Environment1 + (1|group) , 
                     family = "poisson", data = testData)
sr <- simulateResiduals(rightModel, re.form = NA, plot = TRUE)

lme4::glmer() Conditional Model

Use re.form argument in glmer() to resimulate the data.

res2 <- simulateResiduals(fittedModel = rightModel, re.form = NULL, plot = TRUE)

glmmTMB Marginal Models

The implementation is a very recent addition and lacks clear documention (see the entire discussion).

data("sleepstudy", package = "lme4")
m1 <- glmmTMB(Reaction ~ Days  + (Days|Subject), sleepstudy)
m1_res <- simulateResiduals(m1, plot = TRUE)

Marginal Model

Set the simulation conditions for a glmmTMB object. The object is altered in place, hence why it’s copied here.

m2 <- m1
set_simcodes(m2$obj, val = "fix")
m2_res <- simulateResiduals(m1, plot = TRUE)

Other Functionality

(Implemented in DHARMa)

• Uniformity test
• Quantile location tesy
• Outlier test
• Tests for spatial and temporal autocorrelation
• Test for zero-inflation
• Over/Under dispersion tests

Best Practices

• Feng et al (2018) demonstrated for many simulated scenarios that quantile residuals are a better diagnostic tool for checking model distribution, model parameters, over/under dispersion, and overall lack of fit of GLMMs than ‘standard’ residuals.

• The plot of residuals versus fitted is the most reliable indicator of a properly specified model

• We do not know their overall sensitivity of this time to modest model misspecifications

• Results from quantiles tests for uniformity should be treated with caution - no pattern proves a model is correct; likewise, nonrandom matterns are not always a deal-breaker

• Outliers identified are warnings of possible outliers and lack quantitative information

Final Thoughts

• There is a lack of information on using quantile residuals in a mixed model context and the uses of conditional versus marginal residuals
• RTFM: read the DHARMa vignette and when needed, function documentation. Not all all R code documentation is well written and helpful, but this package is very helpful.
• The quantile residuals described herein do not apply to multinomial models. An extension of quantile residuals for these models has been developed, but it lacks an easy implememtation.

Binomial Example

\[y_i \sim binomial(n, p_i) \]

\[ \text{logit} (p_i) = \beta_0 + \beta_1x_i \] \[r_i = \Phi^{-1}(u_i)\]

\[ u_i \sim Uniform(a_i, b_j] \]

\[a_i = lim_{y \uparrow y_i} F(y; n, \hat{p_i}) =\displaystyle\sum_{i=0}^{\lfloor k-1 \rfloor} \begin{pmatrix} n \\ i \end{pmatrix} p_i(1-p)^{n-i}\]

\[b_i = F(y_i; n, \hat{p_i}) =\displaystyle\sum_{i=0}^{\lfloor k \rfloor} \begin{pmatrix} n \\ i \end{pmatrix} p_i(1-p)^{n-i}\]

Consider Groups

This is the ‘wrong model’ (missing ‘environment’), but the plot looks okay

wrongModel2 <- glmer(observedResponse ~ 1 + (1|group), 
                     family = "poisson", data = testData)
res2 <- simulateResiduals(wrongModel2, plot = TRUE)

Consider Groups

When plotted by the missing group:

plotResiduals(res2, form = testData$Environment1)

Standard Residual Plots